#include "spaceDiscr.h"
#include <iostream>
using namespace std;

SpaceDiscr::SpaceDiscr(const Config &cfg, const Geometry &geo,
                       Fluid &flow)
    : si{geo.si}, sj{geo.sj}, vol{geo.vol}, dvE{flow.GetdvEuler()},
      dvNS{flow.GetdvNS()}, cv{flow.GetCV()}, press{flow.GetPressure()},
      BoundDescrp{geo.BoundDescrp}, limCoef{cfg.GetLimiterCoef()},
      epsEntr{cfg.GetEntropyFixCoef()},
      convFluxType{cfg.GetConvectFluxType()}, equsType{cfg.GetEquationType()},
      pExtrNum{cfg.GetPressureExtrNum()}, imax{geo.imax}, jmax{geo.jmax},
      il{geo.il}, jl{geo.jl}, i2{geo.i2}, j2{geo.j2}

{
  int imaxp1{imax + 1}, jmaxp1{jmax + 1};
  diss.ReConstruct(imaxp1, jmaxp1);
  rhs.ReConstruct(imaxp1, jmaxp1);

  if (equsType == Equation_t::NavierStokes)
  {
    gradFi.ReConstruct(imaxp1, jmaxp1);
    gradFj.ReConstruct(imaxp1, jmaxp1);
    u.ReConstruct(imaxp1, jmaxp1);
    v.ReConstruct(imaxp1, jmaxp1);
  }
}

void SpaceDiscr::GradientInitial(void)
{
  // 初始化
  for (int j = 0; j <= jmax; ++j)
  {
    for (int i = 0; i <= imax; ++i)
    {
      const ConsVar_t &tmpCv{cv(i, j)};
      u(i, j) = tmpCv.rhou / tmpCv.rho;
      v(i, j) = tmpCv.rhov / tmpCv.rho;
    }
  }
  gradFi = Grad_t{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  gradFj = Grad_t{0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
}

void SpaceDiscr::GradientFaceI(void)
{
  real_t sx, sy, uav, vav, Tav, rvol;
  real_t fgx[3], fgy[3];

  for (int j = 2; j <= jl; ++j)
  {
    for (int i = 2; i <= i2; ++i)
    {
      // 辅助控制体的右界面
      sx = 0.5 * (si(i, j).sx + si(i + 1, j).sx);
      sy = 0.5 * (si(i, j).sy + si(i + 1, j).sy);
      fgx[0] = u(i, j) * sx;
      fgy[0] = u(i, j) * sy;
      fgx[1] = v(i, j) * sx;
      fgy[1] = v(i, j) * sy;
      fgx[2] = dvNS(i, j).T * sx;
      fgy[2] = dvNS(i, j).T * sy;
      // 向外为正
      Grad_t &grIJ{gradFi(i, j)};
      grIJ.ux -= fgx[0];
      grIJ.uy -= fgy[0];
      grIJ.vx -= fgx[1];
      grIJ.vy -= fgy[1];
      grIJ.Tx -= fgx[2];
      grIJ.Ty -= fgy[2];

      Grad_t &grIp1J{gradFi(i + 1, j)};
      grIp1J.ux += fgx[0];
      grIp1J.uy += fgy[0];
      grIp1J.vx += fgx[1];
      grIp1J.vy += fgy[1];
      grIp1J.Tx += fgx[2];
      grIp1J.Ty += fgy[2];

      // 底部的面  i = 2 的辅助控制体 底部和上部面只计算了一半
      sx = 0.5 * (sj(i, j).sx + sj(i - 1, j).sx);
      sy = 0.5 * (sj(i, j).sy + sj(i - 1, j).sy);

      uav = 0.25 * (u(i, j) + u(i - 1, j) + u(i, j - 1) + u(i - 1, j - 1));
      vav = 0.25 * (v(i, j) + v(i - 1, j) + v(i, j - 1) + v(i - 1, j - 1));
      Tav = 0.25 *
            (dvNS(i, j).T + dvNS(i - 1, j).T + dvNS(i, j - 1).T + dvNS(i - 1, j - 1).T);
      fgx[0] = uav * sx;
      fgy[0] = uav * sy;
      fgx[1] = vav * sx;
      fgy[1] = vav * sy;
      fgx[2] = Tav * sx;
      fgy[2] = Tav * sy;

      grIJ.ux += fgx[0];
      grIJ.uy += fgy[0];
      grIJ.vx += fgx[1];
      grIJ.vy += fgy[1];
      grIJ.Tx += fgx[2];
      grIJ.Ty += fgy[2];

      Grad_t &grIJm1{gradFi(i, j - 1)};
      grIJm1.ux -= fgx[0];
      grIJm1.uy -= fgy[0];
      grIJm1.vx -= fgx[1];
      grIJm1.vy -= fgy[1];
      grIJm1.Tx -= fgx[2];
      grIJm1.Ty -= fgy[2];
    }

    //  i= 2 边界
    sx = si(2, j).sx;
    sy = si(2, j).sy;
    fgx[0] = 0.5 * (u(1, j) + u(2, j)) * sx;
    fgy[0] = 0.5 * (u(1, j) + u(2, j)) * sy;
    fgx[1] = 0.5 * (v(1, j) + v(2, j)) * sx;
    fgy[1] = 0.5 * (v(1, j) + v(2, j)) * sy;
    fgx[2] = 0.5 * (dvNS(1, j).T + dvNS(2, j).T) * sx;
    fgy[2] = 0.5 * (dvNS(1, j).T + dvNS(2, j).T) * sy;

    Grad_t &gr2J{gradFi(2, j)};
    gr2J.ux += fgx[0];
    gr2J.uy += fgy[0];
    gr2J.vx += fgx[1];
    gr2J.vy += fgy[1];
    gr2J.Tx += fgx[2];
    gr2J.Ty += fgy[2];

    //  i = il;
    sx = si(il, j).sx;
    sy = si(il, j).sy;
    fgx[0] = 0.5 * (u(i2, j) + u(il, j)) * sx;
    fgy[0] = 0.5 * (u(i2, j) + u(il, j)) * sy;
    fgx[1] = 0.5 * (v(i2, j) + v(il, j)) * sx;
    fgy[1] = 0.5 * (v(i2, j) + v(il, j)) * sy;
    fgx[2] = 0.5 * (dvNS(i2, j).T + dvNS(il, j).T) * sx;
    fgy[2] = 0.5 * (dvNS(i2, j).T + dvNS(il, j).T) * sy;

    Grad_t &grIlJ{gradFi(il, j)};
    grIlJ.ux -= fgx[0];
    grIlJ.uy -= fgy[0];
    grIlJ.vx -= fgx[1];
    grIlJ.vy -= fgy[1];
    grIlJ.Tx -= fgx[2];
    grIlJ.Ty -= fgy[2];

    // i= il 的 底部面  大概认为两侧相同
    sx = 0.5 * sj(i2, j).sx;
    sy = 0.5 * sj(i2, j).sy;
    uav = 0.25 * (u(i2, j) + u(il, j) + u(i2, j - 1) + u(il, j - 1));
    vav = 0.25 * (v(i2, j) + v(il, j) + v(i2, j - 1) + v(il, j - 1));
    Tav =
        0.25 * (dvNS(i2, j).T + dvNS(il, j).T + dvNS(i2, j - 1).T + dvNS(il, j - 1).T);

    fgx[0] = uav * sx;
    fgy[0] = uav * sy;
    fgx[1] = vav * sx;
    fgy[1] = vav * sy;
    fgx[2] = Tav * sx;
    fgy[2] = Tav * sy;

    grIlJ.ux += fgx[0];
    grIlJ.uy += fgy[0];
    grIlJ.vx += fgx[1];
    grIlJ.vy += fgy[1];
    grIlJ.Tx += fgx[2];
    grIlJ.Ty += fgy[2];

    Grad_t &grIlJm1{gradFi(il, j - 1)};
    grIlJm1.ux -= fgx[0];
    grIlJm1.uy -= fgy[0];
    grIlJm1.vx -= fgx[1];
    grIlJm1.vy -= fgy[1];
    grIlJm1.Tx -= fgx[2];
    grIlJm1.Ty -= fgy[2];
  }

  // 除以辅助控制体的体积即为面的梯度
  for (int j = 2; j <= j2; ++j)
  {
    // 内部的
    for (int i = 3; i <= i2; ++i)
    {
      rvol = 2.0 / (vol(i, j) + vol(i - 1, j));
      Grad_t &grIJ{gradFi(i, j)};
      grIJ.ux *= rvol;
      grIJ.uy *= rvol;
      grIJ.vx *= rvol;
      grIJ.vy *= rvol;
      grIJ.Tx *= rvol;
      grIJ.Ty *= rvol;
    }
    // 边界上的体积为 0.5vol
    rvol = 2.0 / vol(2, j);
    Grad_t &gr2J{gradFi(2, j)};
    gr2J.ux *= rvol;
    gr2J.uy *= rvol;
    gr2J.vx *= rvol;
    gr2J.vy *= rvol;
    gr2J.Tx *= rvol;
    gr2J.Ty *= rvol;

    rvol = 2.0 / vol(i2, j);
    Grad_t &grIlJ{gradFi(il, j)};
    grIlJ.ux *= rvol;
    grIlJ.uy *= rvol;
    grIlJ.vx *= rvol;
    grIlJ.vy *= rvol;
    grIlJ.Tx *= rvol;
    grIlJ.Ty *= rvol;
  }
}

void SpaceDiscr::GradientFaceJ(void)
{
  real_t sx, sy, uav, vav, Tav, rvol;
  real_t fgx[3], fgy[3];

  for (int i = 2; i <= il; ++i)
  {
    for (int j = 2; j <= j2; ++j)
    {
      // bottom face of auxiliary control volume
      sx = 0.5 * (sj(i, j).sx + sj(i, j + 1).sx);
      sy = 0.5 * (sj(i, j).sy + sj(i, j + 1).sy);
      fgx[0] = u(i, j) * sx;
      fgy[0] = u(i, j) * sy;
      fgx[1] = v(i, j) * sx;
      fgy[1] = v(i, j) * sy;
      fgx[2] = dvNS(i, j).T * sx;
      fgy[2] = dvNS(i, j).T * sy;
      // 向外为正
      Grad_t &grIJ{gradFj(i, j)};
      grIJ.ux -= fgx[0];
      grIJ.uy -= fgy[0];
      grIJ.vx -= fgx[1];
      grIJ.vy -= fgy[1];
      grIJ.Tx -= fgx[2];
      grIJ.Ty -= fgy[2];

      Grad_t &grIJp1{gradFj(i, j + 1)};
      grIJp1.ux += fgx[0];
      grIJp1.uy += fgy[0];
      grIJp1.vx += fgx[1];
      grIJp1.vy += fgy[1];
      grIJp1.Tx += fgx[2];
      grIJp1.Ty += fgy[2];

      // left face of auxiliary control volume
      sx = 0.5 * (si(i, j).sx + si(i, j - 1).sx);
      sy = 0.5 * (si(i, j).sy + si(i, j - 1).sy);

      uav = 0.25 * (u(i, j) + u(i - 1, j) + u(i, j - 1) + u(i - 1, j - 1));
      vav = 0.25 * (v(i, j) + v(i - 1, j) + v(i, j - 1) + v(i - 1, j - 1));
      Tav = 0.25 *
            (dvNS(i, j).T + dvNS(i - 1, j).T + dvNS(i, j - 1).T + dvNS(i - 1, j - 1).T);
      fgx[0] = uav * sx;
      fgy[0] = uav * sy;
      fgx[1] = vav * sx;
      fgy[1] = vav * sy;
      fgx[2] = Tav * sx;
      fgy[2] = Tav * sy;

      grIJ.ux += fgx[0];
      grIJ.uy += fgy[0];
      grIJ.vx += fgx[1];
      grIJ.vy += fgy[1];
      grIJ.Tx += fgx[2];
      grIJ.Ty += fgy[2];

      Grad_t &grIm1J{gradFj(i - 1, j)};
      grIm1J.ux -= fgx[0];
      grIm1J.uy -= fgy[0];
      grIm1J.vx -= fgx[1];
      grIm1J.vy -= fgy[1];
      grIm1J.Tx -= fgx[2];
      grIm1J.Ty -= fgy[2];
    }
    //  j= 2 边界
    sx = sj(i, 2).sx;
    sy = sj(i, 2).sy;
    fgx[0] = 0.5 * (u(i, 1) + u(i, 2)) * sx;
    fgy[0] = 0.5 * (u(i, 1) + u(i, 2)) * sy;
    fgx[1] = 0.5 * (v(i, 1) + v(i, 2)) * sx;
    fgy[1] = 0.5 * (v(i, 1) + v(i, 2)) * sy;
    fgx[2] = 0.5 * (dvNS(i, 1).T + dvNS(i, 2).T) * sx;
    fgy[2] = 0.5 * (dvNS(i, 1).T + dvNS(i, 2).T) * sy;

    Grad_t &grI2{gradFj(i, 2)};
    grI2.ux += fgx[0];
    grI2.uy += fgy[0];
    grI2.vx += fgx[1];
    grI2.vy += fgy[1];
    grI2.Tx += fgx[2];
    grI2.Ty += fgy[2];

    //  i = il;
    sx = sj(i, jl).sx;
    sy = sj(i, jl).sy;
    fgx[0] = 0.5 * (u(i, j2) + u(i, jl)) * sx;
    fgy[0] = 0.5 * (u(i, j2) + u(i, jl)) * sy;
    fgx[1] = 0.5 * (v(i, j2) + v(i, jl)) * sx;
    fgy[1] = 0.5 * (v(i, j2) + v(i, jl)) * sy;
    fgx[2] = 0.5 * (dvNS(i, j2).T + dvNS(i, jl).T) * sx;
    fgy[2] = 0.5 * (dvNS(i, j2).T + dvNS(i, jl).T) * sy;

    Grad_t &grIJl{gradFj(i, jl)};
    grIJl.ux -= fgx[0];
    grIJl.uy -= fgy[0];
    grIJl.vx -= fgx[1];
    grIJl.vy -= fgy[1];
    grIJl.Tx -= fgx[2];
    grIJl.Ty -= fgy[2];

    // i= il 的 底部面  大概认为两侧相同
    sx = 0.5 * si(i, j2).sx;
    sy = 0.5 * si(i, j2).sy;
    uav = 0.25 * (u(i, j2) + u(i, jl) + u(i - 1, j2) + u(i - 1, jl));
    vav = 0.25 * (v(i, jl) + v(i, j2) + v(i - 1, j2) + v(i - 1, jl));
    Tav =
        0.25 * (dvNS(i, jl).T + dvNS(i, j2).T + dvNS(i - 1, j2).T + dvNS(i - 1, jl).T);

    fgx[0] = uav * sx;
    fgy[0] = uav * sy;
    fgx[1] = vav * sx;
    fgy[1] = vav * sy;
    fgx[2] = Tav * sx;
    fgy[2] = Tav * sy;

    grIJl.ux += fgx[0];
    grIJl.uy += fgy[0];
    grIJl.vx += fgx[1];
    grIJl.vy += fgy[1];
    grIJl.Tx += fgx[2];
    grIJl.Ty += fgy[2];

    Grad_t &grIm1Jl{gradFj(i - 1, jl)};
    grIm1Jl.ux -= fgx[0];
    grIm1Jl.uy -= fgy[0];
    grIm1Jl.vx -= fgx[1];
    grIm1Jl.vy -= fgy[1];
    grIm1Jl.Tx -= fgx[2];
    grIm1Jl.Ty -= fgy[2];
  }

  // 除以辅助控制体的体积即为面的梯度
  for (int i = 2; i <= i2; ++i)
  {
    // 内部的
    for (int j = 3; j <= j2; ++j)
    {
      rvol = 2.0 / (vol(i, j) + vol(i, j - 1));
      Grad_t &grIJ{gradFj(i, j)};
      grIJ.ux *= rvol;
      grIJ.uy *= rvol;
      grIJ.vx *= rvol;
      grIJ.vy *= rvol;
      grIJ.Tx *= rvol;
      grIJ.Ty *= rvol;
    }
    // 边界上的体积为 0.5vol
    rvol = 2.0 / vol(i, 2);
    Grad_t &grI2{gradFj(i, 2)};
    grI2.ux *= rvol;
    grI2.uy *= rvol;
    grI2.vx *= rvol;
    grI2.vy *= rvol;
    grI2.Tx *= rvol;
    grI2.Ty *= rvol;

    rvol = 2.0 / vol(i, j2);
    Grad_t &grIJl{gradFj(i, jl)};
    grIJl.ux *= rvol;
    grIJl.uy *= rvol;
    grIJl.vx *= rvol;
    grIJl.vy *= rvol;
    grIJl.Tx *= rvol;
    grIJl.Ty *= rvol;
  }
}

void SpaceDiscr::FluxViscous(void)
{
  // 是否要置零  需要置零  但 roe1耗散项也可能置零
  // ================================================
  // 注意： 需要在外部置零 注意与roe耗散项的配合
  constexpr real_t two3{2.0 / 3.0};
  int jj2, iil, m1, im1, jm1, ii, jj;
  real_t uav, vav, muav, kaav, tauxx, tauyy, tauxy, phix, phiy, fv[3];
  const int index2[2]{j2, i2}, indexl[2]{il, jl};
  for (int k = 0; k < 2; ++k)
  {
    // jj2 = (k == 0) ? j2 : i2;
    // iil = (k == 0) ? il : jl;
    jj2 = index2[k];
    iil = indexl[k];
    // const int(&trans)[4]{(k == 0) ? transJ : transI};
    Transij trans(k != 0);
    // 下面的不能使用 if 因为引用必须在定义时绑定
    const VarArray<Area_t, 2, 0> &S{(k == 0) ? si : sj};
    const VarArray<Grad_t, 2, 0> &gradF{(k == 0) ? gradFi : gradFj};

    for (int j = 2; j <= jj2; ++j)
    {
      for (int i = 2; i <= iil; ++i)
      {
        // cout << "j = " << j << "\t" << "i = " << i << endl;
        m1 = i - 1;
        trans(m1, j, im1, jm1);
        // 错了   i j 也要交换
        // cout << "jm1 = " << jm1 << "\t" << "im1 = " << im1 << endl;
        trans(i, j, ii, jj);
        uav = 0.5 * (u(im1, jm1) + u(ii, jj));
        vav = 0.5 * (v(im1, jm1) + v(ii, jj));
        muav = 0.5 * (dvNS(im1, jm1).mu + dvNS(ii, jj).mu);
        kaav = 0.5 * (dvNS(im1, jm1).ka + dvNS(ii, jj).ka);

        const Grad_t &grd{gradF(ii, jj)};
        tauxx = two3 * muav * (2.0 * grd.ux - grd.vy);
        tauyy = two3 * muav * (2.0 * grd.vy - grd.ux);
        tauxy = muav * (grd.uy + grd.vx);
        phix = uav * tauxx + vav * tauxy + kaav * grd.Tx;
        phiy = uav * tauxy + vav * tauyy + kaav * grd.Ty;

        const Area_t &sij{S(ii, jj)};
        fv[0] = sij.sx * tauxx + sij.sy * tauxy;
        fv[1] = sij.sx * tauxy + sij.sy * tauyy;
        fv[2] = sij.sx * phix + sij.sy * phiy;

        diss(ii, jj).rhou += fv[0];
        diss(ii, jj).rhov += fv[1];
        diss(ii, jj).rhoE += fv[2];
        diss(im1, jm1).rhou -= fv[0];
        diss(im1, jm1).rhov -= fv[1];
        diss(im1, jm1).rhoE -= fv[2];
      }
    }
  }
}

// template<typename T>
//   void DissipRoe1(const VarArray<T, 2, 0>& dv);
// template <typename T>
void SpaceDiscr::DissipRoe1(void)
{
  real_t ds, nx, ny, rrho, rl, ul, vl, pl, hl, rr, ur, vr, pr, hr, rav, dd, dd1,
      uav, vav, hav, q2a, c2a, cav, uv, du, h1, h2, h3, h4, h5, delta, eabs1,
      eabs2, eabs4;
  real_t fd[4];

  int im1, jm1, jj2, iil, m1, ii, jj;
  const int index2[2]{j2, i2}, indexl[2]{il, jl};

  for (int k = 0; k < 2; ++k)
  {
    jj2 = index2[k];
    iil = indexl[k];
    // const int(&trans)[4]{(k == 0) ? transJ : transI};
    Transij trans(k != 0);
    // 下面的不能使用 if 因为引用必须在定义时绑定
    const VarArray<Area_t, 2, 0> &S{(k == 0) ? si : sj};

    // 互换 i j 的角色来减少代码量和重复
    for (int j = 2; j <= jj2; ++j)
    {
      for (int i = 2; i <= iil; ++i)
      {
        m1 = i - 1;
        // i j 的角色也互换了
        trans(i, j, ii, jj);
        const Area_t &sij{S(ii, jj)};
        ds = sqrt(sij.sx * sij.sx + sij.sy * sij.sy);
        nx = sij.sx / ds;
        ny = sij.sy / ds;

        trans(m1, j, im1, jm1);
        // ! --- left and right state (at I-1/2)
        const ConsVar_t &cvIm1Jm1{cv(im1, jm1)};
        rrho = 1.0 / cvIm1Jm1.rho;
        rl = cvIm1Jm1.rho;
        ul = cvIm1Jm1.rhou * rrho;
        vl = cvIm1Jm1.rhov * rrho;
        pl = press(im1, jm1);
        //   总焓
        hl = (pl + cvIm1Jm1.rhoE) * rrho;

        const ConsVar_t &cvIJ{cv(ii, jj)};
        rrho = 1.0 / cvIJ.rho;
        rr = cvIJ.rho;
        ur = cvIJ.rhou * rrho;
        vr = cvIJ.rhov * rrho;
        pr = press(ii, jj);
        hr = (pr + cvIJ.rhoE) * rrho;

        // ! --- Roe's average
        rav = sqrt(rl * rr);
        dd = rav / rl;
        dd1 = 1.0 / (1.0 + dd);
        uav = (ul + dd * ur) * dd1;
        vav = (vl + dd * vr) * dd1;
        hav = (hl + dd * hr) * dd1;
        q2a = 0.5 * (uav * uav + vav * vav);
        c2a = gam1 * (hav - q2a);
        cav = sqrt(c2a);
        uv = uav * nx + vav * ny;
        //   ! 面两侧速度矢量的差值在面法向方向上的投影;
        du = (ur - ul) * nx + (vr - vl) * ny;

        // ! --- eigenvalues
        h1 = fabs(uv - cav);
        h2 = fabs(uv);
        h4 = fabs(uv + cav);
        // epsEntr 为熵修正系数
        // 书上是设为当地音速的 1/10  1/20  但实现的时候是 h4
        // 在这里我改为了 cav  h4 不太好 有可能很小
        delta = epsEntr * h4;
        EntropyCorrect(h1, delta, eabs1);
        EntropyCorrect(h2, delta, eabs2);
        EntropyCorrect(h4, delta, eabs4);

        // ! --- upwind fluxes

        h1 = rav * cav * du;
        h2 = eabs1 * (pr - pl - h1) / (2.0 * c2a);
        h3 = eabs2 * (rr - rl - (pr - pl) / c2a);
        h4 = eabs2 * rav;
        h5 = eabs4 * (pr - pl + h1) / (2.0 * c2a);

        fd[0] = h2 + h3 + h5;
        fd[1] = h2 * (uav - cav * nx) + h3 * uav + h4 * (ur - ul - du * nx) +
                h5 * (uav + cav * nx);
        fd[2] = h2 * (vav - cav * ny) + h3 * vav + h4 * (vr - vl - du * ny) +
                h5 * (vav + cav * ny);
        fd[3] = h2 * (hav - cav * uv) + h3 * q2a +
                h4 * (uav * (ur - ul) + vav * (vr - vl) - uv * du) +
                h5 * (hav + cav * uv);

        //             ! --- dissipation term
        //   ! 与 dissipcentral 是同样的
        //   ! 耗散项是流进为正  流出为负  与 flux 异号
        // 不要忘了 0.5   血的教训
        ds *= 0.5;
        diss(ii, jj).rho -= fd[0] * ds;
        diss(ii, jj).rhou -= fd[1] * ds;
        diss(ii, jj).rhov -= fd[2] * ds;
        diss(ii, jj).rhoE -= fd[3] * ds;

        diss(im1, jm1).rho += fd[0] * ds;
        diss(im1, jm1).rhou += fd[1] * ds;
        diss(im1, jm1).rhov += fd[2] * ds;
        diss(im1, jm1).rhoE += fd[3] * ds;
      }
    }
  }

  // cout << " DissipRoe1 diss(4,4) = " << diss(4, 4).rho << "\n"
  //      << diss(4, 4).rhou << "\n"
  //      << diss(4, 4).rhov << "\n"
  //      << diss(4, 4).rhoE << "\n";
}
// 没用计算边界面的通量
// template<typename T>
//   void FluxRoe1(const VarArray<T, 2, 0>& dv);
// template <typename T>
void SpaceDiscr::FluxRoe1(void)
{
  for (int j = 2; j <= j2; ++j)
  {
    for (int i = 2; i <= i2; ++i)
    {
      rhs(i, j).rho = -diss(i, j).rho;
      rhs(i, j).rhou = -diss(i, j).rhou;
      rhs(i, j).rhov = -diss(i, j).rhov;
      rhs(i, j).rhoE = -diss(i, j).rhoE;
    }
  }

  real_t rhl, rhr, qsl, qsr, pav, fc[4];
  for (int j = 2; j <= j2; ++j)
  {
    // flux in i-direction (except through boundary)
    for (int i = 3; i <= i2; ++i)
    {
      // ! rho * H
      const ConsVar_t &cvIJ{cv(i, j)}, &cvIm1J{cv(i - 1, j)};
      const Area_t &sij{si(i, j)};
      rhl = press(i - 1, j) + cvIm1J.rhoE;
      //   ! V .* S   体积流量;
      qsl = (cvIm1J.rhou * sij.sx + cvIm1J.rhov * sij.sy) / cvIm1J.rho;
      rhr = press(i, j) + cvIJ.rhoE;
      qsr = (cvIJ.rhou * sij.sx + cvIJ.rhov * sij.sy) / cvIJ.rho;

      pav = 0.5 * (press(i - 1, j) + press(i, j));
      fc[0] = 0.5 * (qsl * cvIm1J.rho + qsr * cvIJ.rho);
      fc[1] = 0.5 * (qsl * cvIm1J.rhou + qsr * cvIJ.rhou) + sij.sx * pav;
      fc[2] = 0.5 * (qsl * cvIm1J.rhov + qsr * cvIJ.rhov) + sij.sy * pav;
      fc[3] = 0.5 * (qsl * rhl + qsr * rhr);

      rhs(i, j).rho += fc[0];
      rhs(i, j).rhou += fc[1];
      rhs(i, j).rhov += fc[2];
      rhs(i, j).rhoE += fc[3];

      rhs(i - 1, j).rho -= fc[0];
      rhs(i - 1, j).rhou -= fc[1];
      rhs(i - 1, j).rhov -= fc[2];
      rhs(i - 1, j).rhoE -= fc[3];
    }

    // flux in j-direction (except through boundary)
    if (j > 2)
    {
      // 有预测机制  不用担心性能问题  瞎考虑那么多，都没用
      // 还是这个强  i 一直在内循环里
      for (int i = 2; i <= i2; ++i)
      {
        const ConsVar_t &cvIJ{cv(i, j)}, &cvIJm1{cv(i, j - 1)};
        const Area_t &sij{sj(i, j)};
        rhl = press(i, j - 1) + cvIJm1.rhoE;
        qsl = (cvIJm1.rhou * sij.sx + cvIJm1.rhov * sij.sy) / cvIJm1.rho;
        rhr = press(i, j) + cvIJ.rhoE;
        qsr = (cvIJ.rhou * sij.sx + cvIJ.rhov * sij.sy) / cvIJ.rho;

        pav = 0.5 * (press(i, j - 1) + press(i, j));
        fc[0] = 0.5 * (qsl * cvIJm1.rho + qsr * cvIJ.rho);
        fc[1] = 0.5 * (qsl * cvIJm1.rhou + qsr * cvIJ.rhou) + sij.sx * pav;
        fc[2] = 0.5 * (qsl * cvIJm1.rhov + qsr * cvIJ.rhov) + sij.sy * pav;
        fc[3] = 0.5 * (qsl * rhl + qsr * rhr);

        rhs(i, j).rho += fc[0];
        rhs(i, j).rhou += fc[1];
        rhs(i, j).rhov += fc[2];
        rhs(i, j).rhoE += fc[3];

        rhs(i, j - 1).rho -= fc[0];
        rhs(i, j - 1).rhou -= fc[1];
        rhs(i, j - 1).rhov -= fc[2];
        rhs(i, j - 1).rhoE -= fc[3];
      }
    }
  }

  // 处理边界
  int nbound = BoundDescrp.size();
  for (int i = 0; i < nbound; ++i)
  {
    switch (BoundDescrp[i].BType)
    {
    case Bound_t::Wall:
    case Bound_t::Symmetry:
      FluxBoundary(BoundDescrp[i]);
      break;

    case Bound_t::Inflow:
    case Bound_t::Outflow:
    case Bound_t::Farfield:
    case Bound_t::InnerBound:
      FluxRoeBoundary(BoundDescrp[i]);
      break;
    default:
      cerr << "do not support noboundary conditions" << endl;
      exit(EXIT_FAILURE);
      break;
    }
  }
}
// wall 和 symmetry
// template <typename T>
void SpaceDiscr::FluxBoundary(const BoundDescrp_t &Bdescrp)
{
  const bool isIface{Bdescrp.isIface};
  // const int(&trans)[4] = isIface ? transI : transJ;
  Transij trans(isIface);
  const int dum1{Bdescrp.dum1}, ins1{Bdescrp.ins1},
      ins2{Bdescrp.ins2};
  int beg{Bdescrp.beg}, end{Bdescrp.end};

  const VarArray<Area_t, 2, 0> &S{isIface ? si : sj};

  int indexBface{(Bdescrp.Bn < 2) ? ins1 : dum1},
      ins3{(Bdescrp.Bn < 2) ? ins2 + 1 : ins2 - 1};
  // 使面的指向朝外
  real_t sgn{(Bdescrp.Bn < 2) ? 1.0 : -1.0};
  real_t sx, sy, pa;
  int id1, jd1, id2, jd2, id3, jd3;
  switch (Bdescrp.BType)
  {
  case Bound_t::Wall:
    if (equsType == Equation_t::NavierStokes)
    {
      for (int i = beg; i <= end; ++i)
      {
        trans(i, indexBface, id1, jd1);
        // ! 保证面法向向外
        sx = sgn * S(id1, jd1).sx;
        sy = sgn * S(id1, jd1).sy;

        trans(i, ins1, id1, jd1);
        pa = press(id1, jd1);
        // 流进还是流出看面的方向
        rhs(id1, jd1).rhou += pa * sx;
        rhs(id1, jd1).rhov += pa * sy;
      }
    }
    else
    {
      for (int i = beg; i <= end; ++i)
      {
        trans(i, indexBface, id1, jd1);
        // ! 保证面法向向外
        sx = sgn * S(id1, jd1).sx;
        sy = sgn * S(id1, jd1).sy;

        trans(i, ins1, id1, jd1);
        trans(i, ins2, id2, jd2);
        if (pExtrNum == 2)
        {
          pa = 0.5 * (3.0 * press(id1, jd1) - press(id2, jd2));
        }
        else
        {
          trans(i, ins3, id3, jd3);
          pa = (15.0 * press(id1, jd1) - 10.0 * press(id2, jd2) +
                3.0 * press(id3, jd3)) /
               8.0;
        }
        // 流进还是流出看面的方向
        rhs(id1, jd1).rhou += pa * sx;
        rhs(id1, jd1).rhov += pa * sy;
      }
    }
    break;

  case Bound_t::Symmetry:
    for (int i = beg; i <= end; ++i)
    {
      trans(i, indexBface, id1, jd1);
      // ! 保证面法向向外
      sx = sgn * S(id1, jd1).sx;
      sy = sgn * S(id1, jd1).sy;

      trans(i, ins1, id1, jd1);
      trans(i, dum1, id2, jd2);
      // (p(ins1) + p(dum1))*0.5
      pa = 0.5 * (press(id1, jd1) + press(id2, jd2));
      // 流进还是流出看面的方向
      rhs(id1, jd1).rhou += pa * sx;
      rhs(id1, jd1).rhov += pa * sy;
    }
    break;
  default:
    cerr << "error: this function only support wall and symmetry" << endl;
    exit(EXIT_FAILURE);
    break;
  }
}
// inflow outflow farfield 周期  cut      这些都是黎曼问题
// 暂时没有考虑二阶精度
// template <typename T>
void SpaceDiscr::FluxRoeBoundary(const BoundDescrp_t &Bdescrp)
{
  const bool isIface{Bdescrp.isIface};
  // const int(&trans)[4] = isIface ? transI : transJ;
  Transij trans(isIface);
  const int dum1{Bdescrp.dum1}, ins1{Bdescrp.ins1}, ins2{Bdescrp.ins2};
  int beg{Bdescrp.beg}, end{Bdescrp.end};

  const VarArray<Area_t, 2, 0> &S{isIface ? si : sj};

  int indexBface{(Bdescrp.Bn < 2) ? ins1 : dum1};
  real_t sgn{(Bdescrp.Bn < 2) ? 1.0 : -1.0};

  real_t sx, sy, rrho, rl, ul, vl, pl, hl, rr, ur, vr, pr, hr, qsrl, qsrr, pav;
  real_t fc[4];

  int idum1, jdum1, iins1, jins1;
  for (int i = beg; i <= end; ++i)
  {
    // idum1 下标 临时借给 face 用
    trans(i, indexBface, idum1, jdum1);
    sx = sgn * S(idum1, jdum1).sx;
    sy = sgn * S(idum1, jdum1).sy;

    if (convFluxType == ConvectFlux_t::Roe1)
    {
      // 利用左右各单元状态确定边界面左右状态
      trans(i, dum1, idum1, jdum1);
      const ConsVar_t &cvIJdum1{cv(idum1, jdum1)};
      rrho = 1.0 / cvIJdum1.rho;
      rl = cvIJdum1.rho;
      ul = cvIJdum1.rhou * rrho;
      vl = cvIJdum1.rhov * rrho;
      pl = press(idum1, jdum1);
      hl = (pl + cvIJdum1.rhoE) * rrho;
      // ! 这里 左和右是根据面的方向判定的
      trans(i, ins1, iins1, jins1);
      const ConsVar_t &cvIJins1{cv(iins1, jins1)};
      rrho = 1.0 / cvIJins1.rho;
      rr = cvIJins1.rho;
      ur = cvIJins1.rhou * rrho;
      vr = cvIJins1.rhov * rrho;
      pr = press(iins1, jins1);
      hr = (pr + cvIJins1.rhoE) * rrho;
    }
    else
    {
      cerr << "error: this function only support Roe1 now" << endl;
      exit(EXIT_FAILURE);
    }

    qsrl = (ul * sx + vl * sy) * rl;
    qsrr = (ur * sx + vr * sy) * rr;
    pav = 0.5 * (pl + pr);
    fc[0] = 0.5 * (qsrl + qsrr);
    fc[1] = 0.5 * (qsrl * ul + qsrr * ur) + sx * pav;
    fc[2] = 0.5 * (qsrl * vl + qsrr * vr) + sy * pav;
    fc[3] = 0.5 * (qsrl * hl + qsrr * hr);

    rhs(iins1, jins1).rho += fc[0];
    rhs(iins1, jins1).rhou += fc[1];
    rhs(iins1, jins1).rhov += fc[2];
    rhs(iins1, jins1).rhoE += fc[3];
  }
}
